{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table border=\"0\">\n",
    "    <tr>\n",
    "        <td>\n",
    "            <img src=\"https://ictd2016.files.wordpress.com/2016/04/microsoft-research-logo-copy.jpg\" style=\"width 30px;\" />\n",
    "             </td>\n",
    "        <td>\n",
    "            <img src=\"https://www.microsoft.com/en-us/research/wp-content/uploads/2016/12/MSR-ALICE-HeaderGraphic-1920x720_1-800x550.jpg\" style=\"width 100px;\"/></td>\n",
    "        </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dynamic Double Machine Learning: Use Cases and Examples\n",
    "\n",
    "Dynamic DoubleML is an extension of the Double ML approach for treatments assigned sequentially over time periods. This estimator will account for treatments that can have causal effects on future outcomes. For more details, see [this paper](https://arxiv.org/abs/2002.07285) or the [EconML docummentation](https://www.pywhy.org/EconML/).\n",
    "\n",
    "For example, the Dynamic DoubleML could be useful in estimating the following causal effects:\n",
    "* the effect of investments on revenue at companies that receive investments at regular intervals ([see more](https://arxiv.org/abs/2103.08390))\n",
    "* the effect of prices on demand in stores where prices of goods change over time\n",
    "* the effect of income on health outcomes in people who receive yearly income\n",
    "\n",
    "The preferred data format is balanced panel data. Each panel corresponds to one entity (e.g. company, store or person) and the different rows in a panel correspond to different time points. Example:\n",
    "\n",
    "||Company|Year|Features|Investment|Revenue|\n",
    "|---|---|---|---|---|---|\n",
    "|1|A|2018|...|\\$1,000|\\$10,000|\n",
    "|2|A|2019|...|\\$2,000|\\$12,000|\n",
    "|3|A|2020|...|\\$3,000|\\$15,000|\n",
    "|4|B|2018|...|\\$0|\\$5,000|\n",
    "|5|B|2019|...|\\$100|\\$10,000|\n",
    "|6|B|2020|...|\\$1,200|\\$7,000|\n",
    "|7|C|2018|...|\\$1,000|\\$20,000|\n",
    "|8|C|2019|...|\\$1,500|\\$25,000|\n",
    "|9|C|2020|...|\\$500|\\$15,000|\n",
    "\n",
    "(Note: when passing the data to the DynamicDML estimator, the \"Company\" column above corresponds to the `groups` argument at fit time. The \"Year\" column above should not be passed in as it will be inferred from the \"Company\" column)\n",
    "\n",
    "If group memebers do not appear together, it is assumed that the first instance of a group in the dataset corresponds to the first period of that group, the second instance of the group corresponds to the second period, etc. Example:\n",
    "\n",
    "||Company|Features|Investment|Revenue|\n",
    "|---|---|---|---|---|\n",
    "|1|A|...|\\$1,000|\\$10,000|\n",
    "|2|B|...|\\$0|\\$5,000\n",
    "|3|C|...|\\$1,000|\\$20,000|\n",
    "|4|A|...|\\$2,000|\\$12,000|\n",
    "|5|B|...|\\$100|\\$10,000|\n",
    "|6|C|...|\\$1,500|\\$25,000|\n",
    "|7|A|...|\\$3,000|\\$15,000|\n",
    "|8|B|...|\\$1,200|\\$7,000|\n",
    "|9|C|...|\\$500|\\$15,000|\n",
    "\n",
    "In this dataset, 1<sup>st</sup> row corresponds to the first period of group `A`, 4<sup>th</sup> row corresponds to the second period of group `A`, etc.\n",
    "\n",
    "In this notebook, we show the performance of the DynamicDML on synthetic and observational data. \n",
    "\n",
    "## Notebook Contents\n",
    "\n",
    "1. [Example Usage with Average Treatment Effects](#1.-Example-Usage-with-Average-Treatment-Effects)\n",
    "2. [Example Usage with Heterogeneous Treatment Effects](#2.-Example-Usage-with-Heterogeneous-Treatment-Effects)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main imports\n",
    "from econml.panel.dml import DynamicDML\n",
    "from econml.tests.dgp import DynamicPanelDGP, add_vlines\n",
    "\n",
    "# Helper imports\n",
    "import numpy as np\n",
    "from sklearn.linear_model import LassoCV, MultiTaskLassoCV\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Example Usage with Average Treatment Effects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 DGP\n",
    "\n",
    "We consider a data generating process from a markovian treatment model. \n",
    "\n",
    "In the example bellow, $T_t\\rightarrow$ treatment(s) at time $t$, $Y_t\\rightarrow$outcome at time $t$, $X_t\\rightarrow$ features and controls at time $t$ (the coefficients $e, f$ will pick the features and the controls).\n",
    "\\begin{align}\n",
    "    X_t =& (\\pi'X_{t-1} + 1) \\cdot A\\, T_{t-1} + B X_{t-1} + \\epsilon_t\\\\\n",
    "    T_t =& \\gamma\\, T_{t-1} + (1-\\gamma) \\cdot D X_t + \\zeta_t\\\\\n",
    "    Y_t =& (\\sigma' X_{t} + 1) \\cdot e\\, T_{t} + f X_t + \\eta_t\n",
    "\\end{align}\n",
    "\n",
    "with $X_0, T_0 = 0$ and $\\epsilon_t, \\zeta_t, \\eta_t \\sim N(0, \\sigma^2)$. Moreover, $X_t \\in R^{n_x}$, $B[:, 0:s_x] \\neq 0$ and $B[:, s_x:-1] = 0$, $\\gamma\\in [0, 1]$, $D[:, 0:s_x] \\neq 0$, $D[:, s_x:-1]=0$, $f[0:s_x]\\neq 0$, $f[s_x:-1]=0$. We draw a single time series of samples of length $n\\_panels \\cdot n\\_periods$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define DGP parameters\n",
    "np.random.seed(123)\n",
    "n_panels = 5000 # number of panels\n",
    "n_periods = 3 # number of time periods in each panel\n",
    "n_treatments = 2 # number of treatments in each period\n",
    "n_x = 100 # number of features + controls\n",
    "s_x = 10 # number of controls (endogeneous variables)\n",
    "s_t = 10 # treatment support size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate data\n",
    "dgp = DynamicPanelDGP(n_periods, n_treatments, n_x).create_instance(\n",
    "            s_x, random_seed=12345)\n",
    "Y, T, X, W, groups = dgp.observational_data(n_panels, s_t=s_t, random_seed=12345)\n",
    "true_effect = dgp.true_effect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 Train Estimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "est = DynamicDML(\n",
    "    model_y=LassoCV(cv=3, max_iter=1000),\n",
    "    model_t=MultiTaskLassoCV(cv=3, max_iter=1000),\n",
    "    cv=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<econml.panel.dml._dml.DynamicDML at 0x19d2abd6a00>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est.fit(Y, T, X=None, W=W, groups=groups)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average effect of default policy: 1.40\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "A scalar was specified but there are multiple treatments; the same value will be used for each treatment.  Consider specifyingall treatments, or using the const_marginal_effect method.\n"
     ]
    }
   ],
   "source": [
    "# Average treatment effect of all periods on last period for unit treatments\n",
    "print(f\"Average effect of default policy: {est.ate():0.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Effect of target policy over baseline policy: 1.40\n"
     ]
    }
   ],
   "source": [
    "# Effect of target policy over baseline policy\n",
    "# Must specify a treatment for each period\n",
    "baseline_policy = np.zeros((1, n_periods * n_treatments))\n",
    "target_policy = np.ones((1, n_periods * n_treatments))\n",
    "eff = est.effect(T0=baseline_policy, T1=target_policy)\n",
    "print(f\"Effect of target policy over baseline policy: {eff[0]:0.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Marginal effect of a treatments in period 1 on period 3 outcome: [0.04000235 0.0701606 ]\n",
      "Marginal effect of a treatments in period 2 on period 3 outcome: [0.31611764 0.23714736]\n",
      "Marginal effect of a treatments in period 3 on period 3 outcome: [0.13108411 0.60656886]\n"
     ]
    }
   ],
   "source": [
    "# Period treatment effects + interpretation\n",
    "for i, theta in enumerate(est.intercept_.reshape(-1, n_treatments)):\n",
    "    print(f\"Marginal effect of a treatments in period {i+1} on period {n_periods} outcome: {theta}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient Results:  X is None, please call intercept_inference to learn the constant!\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>CATE Intercept Results</caption>\n",
       "<tr>\n",
       "             <td></td>             <th>point_estimate</th> <th>stderr</th>  <th>zstat</th>  <th>pvalue</th> <th>ci_lower</th> <th>ci_upper</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T0)$_0$</th>      <td>0.04</td>       <td>0.041</td>  <td>0.977</td>   <td>0.328</td>  <td>-0.027</td>    <td>0.107</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T1)$_0$</th>      <td>0.07</td>       <td>0.04</td>   <td>1.74</td>    <td>0.082</td>   <td>0.004</td>    <td>0.136</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T0)$_1$</th>      <td>0.316</td>      <td>0.036</td>  <td>8.848</td>    <td>0.0</td>    <td>0.257</td>    <td>0.375</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T1)$_1$</th>      <td>0.237</td>      <td>0.036</td>  <td>6.608</td>    <td>0.0</td>    <td>0.178</td>    <td>0.296</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T0)$_2$</th>      <td>0.131</td>      <td>0.003</td> <td>45.665</td>    <td>0.0</td>    <td>0.126</td>    <td>0.136</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T1)$_2$</th>      <td>0.607</td>      <td>0.003</td> <td>210.244</td>   <td>0.0</td>    <td>0.602</td>    <td>0.611</td> \n",
       "</tr>\n",
       "</table><br/><br/><sub>A linear parametric conditional average treatment effect (CATE) model was fitted:<br/>$Y = \\Theta(X)\\cdot T + g(X, W) + \\epsilon$<br/>where for every outcome $i$ and treatment $j$ the CATE $\\Theta_{ij}(X)$ has the form:<br/>$\\Theta_{ij}(X) = \\phi(X)' coef_{ij} + cate\\_intercept_{ij}$<br/>where $\\phi(X)$ is the output of the `featurizer` or $X$ if `featurizer`=None. Coefficient Results table portrays the $coef_{ij}$ parameter vector for each outcome $i$ and treatment $j$. Intercept Results table portrays the $cate\\_intercept_{ij}$ parameter.</sub>"
      ],
      "text/plain": [
       "<class 'econml.utilities.Summary'>\n",
       "\"\"\"\n",
       "                            CATE Intercept Results                            \n",
       "==============================================================================\n",
       "                        point_estimate stderr  zstat  pvalue ci_lower ci_upper\n",
       "------------------------------------------------------------------------------\n",
       "cate_intercept|(T0)$_0$           0.04  0.041   0.977  0.328   -0.027    0.107\n",
       "cate_intercept|(T1)$_0$           0.07   0.04    1.74  0.082    0.004    0.136\n",
       "cate_intercept|(T0)$_1$          0.316  0.036   8.848    0.0    0.257    0.375\n",
       "cate_intercept|(T1)$_1$          0.237  0.036   6.608    0.0    0.178    0.296\n",
       "cate_intercept|(T0)$_2$          0.131  0.003  45.665    0.0    0.126    0.136\n",
       "cate_intercept|(T1)$_2$          0.607  0.003 210.244    0.0    0.602    0.611\n",
       "------------------------------------------------------------------------------\n",
       "\n",
       "<sub>A linear parametric conditional average treatment effect (CATE) model was fitted:\n",
       "$Y = \\Theta(X)\\cdot T + g(X, W) + \\epsilon$\n",
       "where for every outcome $i$ and treatment $j$ the CATE $\\Theta_{ij}(X)$ has the form:\n",
       "$\\Theta_{ij}(X) = \\phi(X)' coef_{ij} + cate\\_intercept_{ij}$\n",
       "where $\\phi(X)$ is the output of the `featurizer` or $X$ if `featurizer`=None. Coefficient Results table portrays the $coef_{ij}$ parameter vector for each outcome $i$ and treatment $j$. Intercept Results table portrays the $cate\\_intercept_{ij}$ parameter.</sub>\n",
       "\"\"\""
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Period treatment effects with confidence intervals\n",
    "est.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "conf_ints = est.intercept__interval(alpha=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 Performance Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3gAAAEzCAYAAABjbqHIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAApUklEQVR4nO3de3TV5Z3v8c/XGEiAXLiJkHAJHWTkThOp1Bs4KNg6x9RedHC1Vo8HOQU9p6elheW1q2Ori47tsHRE2nG0Z4rYUUSLzDBeoFOPtJBIhAJSaARJUAyBhABJSMJz/tgxTcKGBMiPh/3s92st1s7v+T3Z+Qj4XfvL8/x+P3POCQAAAACQ+C7wHQAAAAAA0DVo8AAAAAAgEDR4AAAAABAIGjwAAAAACAQNHgAAAAAEggYPAAAAAAJxoe8Ap6tfv35u2LBhvmMACMihugZJUmZaquckAEJCbQEQleLi4v3Ouf7xziVcgzds2DAVFRX5jgEgICV7qiRJEwZne80BICzUFgBRMbPdJzuXcA0eAHQ1PnwBiAK1BYAPXIMHAAAAAIGgwQOQ9H63o0K/21HhOwaAwFBbAPgQxBbNhoYGlZWVqa6uzncUNEtLS1Nubq5SU7mwHOe/+objviMACBC1BYAPQTR4ZWVlysjI0LBhw2RmvuMkPeecKisrVVZWpry8PN9xAAAAgKQRxBbNuro69e3bl+buPGFm6tu3LyuqAAAAwDkWRIMniebuPMOfBwAAAHDuBdPgna5bnl6nW55e12Xvl5KSogkTJmj06NEaP368Hn/8cR0/7mfvfVFRke69995Tzhk2bJjGjh2rsWPHatSoUbr//vtVX18vSdq1a5fMTA888EDL/P379ys1NVVz586VJD388MP6yU9+Et1/BHAOXZyVpouz0nzHABAYagsAH5K2wetq6enpKikp0ZYtW/T6669r1apV+sEPfuAlS0FBgRYtWtThvDVr1mjz5s1av369SktLNWvWrJZzw4cP18qVK1uO/+3f/k2jR4+OJC/g25icLI3JyfIdA0BgqC0AfEjKBm/FxnJt/LBKf/jggK549C2t2Fjepe9/0UUXacmSJXriiSfknNNVV12lkpKSlvNXXHGFNm3apIcfflh33nmnpkyZouHDh7dpygoLC5Wfn6/Ro0dryZIlLeO9evXS97//feXn52vatGlav359y/e/+uqrkqS1a9fqxhtvlCQdPnxYd9xxh8aOHatx48bppZdeOiFvr169tHjxYq1YsUIHDhyQFGtYL730UhUVFUmSXnjhBX3ta1/r0t8nAAAAAF0r6Rq8FRvLtWD5Zh1rim2fLK+q1YLlm7u8yRs+fLiOHz+uTz75RHfddZeeffZZSdKf/vQn1dfXa9y4cZKk999/X6tXr9b69ev1gx/8QA0NDZKkZ555RsXFxSoqKtKiRYtUWVkpSTpy5IimTJmi4uJiZWRk6P7779frr7+ul19+WQ8++OAJOX74wx8qKytLmzdv1qZNm3TttdfGzZuZmam8vDzt2LGjZezWW2/VsmXLVFZWppSUFA0aNKgrf4uA88aa7Z9ozfZPfMcAEBhqC9A5KzaW64pH31Le/Nf+svhSViyt/I70r1+JvZYV+46ZMCJt8MxshpltN7OdZjb/JHOmmFmJmW0xs99GmUeSFq7ertqGpjZjtQ1NWrh6e5f/LOecJOmrX/2qVq5cqYaGBj3zzDP65je/2TLni1/8orp3765+/frpoosu0r59+yRJixYt0vjx43X55Zdrz549LY1Xt27dNGPGDEnS2LFjdc011yg1NVVjx47Vrl27TsjwxhtvaM6cOS3HvXv37jDvp2bMmKHXX39dzz//vG655ZYz+j0AEkFTk1NTk+t4IgCcBmoL0LFPF1/Kq2rlFFt8+dXyl7Xv1Qek2kopc2Dsdc0jNHmdFNlz8MwsRdKTkq6TVCZpg5m96pzb2mpOtqR/kjTDOfehmV0UVZ5P7a2qPa3xM1VaWqqUlBRddNFFMjNdd911euWVV/TrX/+6ZdujJHXv3r3l65SUFDU2Nmrt2rV64403tG7dOvXo0UNTpkxpeeRAampqyx0qL7jggpbvv+CCC9TY2HhCDudcp+5oWVNTo127dumSSy5RdXW1pFgzmZ+fr3/4h3/Qli1b9Jvf/ObMf0MAAACAduItvvytW6PN+6UBg5uvYU1rfi1ZKuXmn+OEiSfKFbxJknY650qdc8ckLZN0U7s5MyUtd859KEnOucj3MQzKTj+t8TNRUVGh2bNna+7cuS3N1V133aV7771Xl112mfr06XPK76+urlbv3r3Vo0cPvf/++/r9739/xlmuv/56PfHEEy3HBw8ePGHO4cOH9a1vfUuFhYUnrPB95zvf0WOPPaa+ffuecQYAAAAgnniLLIOtQhXHurUd7J4hVe0+R6kSW5QNXo6kPa2Oy5rHWrtEUm8zW2tmxWb2jQjzSJLmTR+p9NSUNmPpqSmaN33kWb1vbW1ty2MSpk2bpuuvv14PPfRQy/n8/HxlZmbqjjvu6PC9ZsyYocbGRo0bN04PPPCALr/88jPOdf/99+vgwYMaM2aMxo8frzVr1rScmzp1qsaMGaNJkyZpyJAhevrpp0/4/tGjR+v222+P+95///d/r9zc3JZfAAAAwOmIt8iyx/VX/27H2g7W10jZQ89RqsRm7a+76rI3NvuqpOnOubuaj78uaZJz7p5Wc56QVCDpbySlS1on6YvOuT+1e69ZkmZJ0pAhQ/J3727bvW/btk2XXnppp7Ot2Fiu7724SceajisnO13zpo9U4cT2vWfX2rt3r6ZMmaL3339fF1yQHPe2Od0/F8CXbR8dkiRdOjDTcxIAIaG2AB379Bq81ts0L0v9QE8MXKUB/S+KrdzV10h1h6Sp97FFs5mZFTvnCuKdi+waPMVW7Aa3Os6VtDfOnP3OuSOSjpjZf0kaL6lNg+ecWyJpiSQVFBScdUdaODFHz6//UJL0wt2Tz/btOvTLX/5S9913nx5//PGkae6ARMKHLwBRoLYAHft0kWXh6u3aW1WrQdnpum36lzSg/+TYNXdVu2Mrd5PvobnrpChX8C5UrFH7G0nlkjZImumc29JqzqWSnpA0XVI3Sesl3eqc++PJ3regoMC1vkmJxErR+Yo/FwAAAKDreVnBc841mtlcSaslpUh6xjm3xcxmN59f7JzbZmb/IWmTpOOSfnGq5g4AovDG1tjjSaaNGuA5CYCQUFsA+BDlFk0551ZJWtVubHG744WSFkaZAwAAAACSAReEAQAAAEAgaPAAAAAAIBA0eF1k3759mjlzpoYPH678/HxNnjxZL7/88jnNsGvXLo0ZMybu+NKlS8/oPX/2s5/p6NGjLce9evU643wAAAAAokWD1wWccyosLNTVV1+t0tJSFRcXa9myZSorKzthbmNj4znPd6oGr6M87Rs8IERD+/bQ0L49fMcAEBhqCwAfIr3JynmrrLjtczUmzDyr52q89dZb6tatm2bPnt0yNnToUN1zT+yZ7s8++6xee+011dXV6ciRI3rxxRd15513qrS0VD169NCSJUs0btw4Pfzww+rVq5e++93vSpLGjBmjlStXSpJuuOEGXXnllXrnnXeUk5OjV155Renp6SouLtadd96pHj166Morr4ybb/78+dq2bZsmTJig22+/Xb17926T58EHH9RPfvKTlp81d+5cFRQU6NChQ9q7d6+mTp2qfv36ac2aNZKk++67TytXrlR6erpeeeUVDRjA3cGQ2EYMyPAdAUCAqC0AfEi+FbyyYmnNI1JtpZQ5MPa65pHY+BnasmWLPvvZz55yzrp16/Tcc8/prbfe0kMPPaSJEydq06ZN+tGPfqRvfOMbHf6MHTt2aM6cOdqyZYuys7P10ksvSZLuuOMOLVq0SOvWrTvp9z766KO66qqrVFJSom9/+9sn5DmZe++9V4MGDdKaNWtamrsjR47o8ssv13vvvaerr75aP//5zzvMDpzvGpuOq7HpuO8YAAJDbQHgQ/I1eCVLpbRMKS1Lsgtir2mZsfEuMmfOHI0fP16XXXZZy9h1112nPn36SJLefvttff3rX5ckXXvttaqsrFR1dfUp3zMvL08TJkyQJOXn52vXrl2qrq5WVVWVrrnmGklqec/OaJ3ndHTr1k033nhjmxxAolu7vUJrt1f4jgEgMNQWAD4kX4NXtVvq3m7LRPeM2PgZGj16tN59992W4yeffFJvvvmmKir+UtR79uzZ8rVz7oT3MDNdeOGFOn78L//SV1dX95eI3bu3fJ2SkqLGxkY552RmZ5S5dZ5T/dz2UlNTW37mpzkAAAAAnB+Sr8HLHirV17Qdq6+JjZ+ha6+9VnV1dXrqqadaxk51Y5Krr75av/rVryRJa9euVb9+/ZSZmalhw4a1NIrvvvuuPvjgg1P+3OzsbGVlZentt9+WpJb3bC8jI0M1NTVxz0mx6wW3bt2q+vp6VVdX68033+z09wIAAAA4fyRfgzdhplR3SKqrltzx2Gvdodj4GTIzrVixQr/97W+Vl5enSZMm6fbbb9djjz0Wd/7DDz+soqIijRs3TvPnz9dzzz0nSfryl7+sAwcOaMKECXrqqad0ySWXdPiz/+Vf/kVz5szR5MmTlZ6eHnfOuHHjdOGFF2r8+PH66U9/esL5wYMH62tf+5rGjRun2267TRMnTmw5N2vWLN1www2aOnVqZ34rAAAAAHhk8bYLns8KCgpcUVFRm7Ft27bp0ksv7fybdPFdNBHfaf+5AJ68sXWfJGnaKO4IC6DrUFsARMXMip1zBfHOJedjEnLzaegAtBjev2fHkwDgNFFbAPiQnA0eALQyvH8v3xEABIjaAsCH5LsGDwDaqWtoUl1Dk+8YAAJDbQHgQzANXqJdSxg6/jyQSN7esV9v79jvOwaAwFBbAPgQRIOXlpamyspKmorzhHNOlZWVSktL8x0FAAAASCpBXIOXm5ursrKyNg8Wh19paWnKzc31HQMAAABIKkE0eKmpqcrLy/MdAwAAAAC8CmKLJgAAAAAgkBU8ADgbIwZwK3MAXY/aAsAHGjwASW9oXx5GDKDrUVsA+MAWTQBJ70h9o47UN/qOASAw1BYAPtDgAUh66/5cqXV/rvQdA0BgqC0AfKDBAwAAAIBA0OABAAAAQCBo8AAAAAAgEDR4AAAAABAIHpMAIOn99cAM3xEABIjaAsAHGjwASS+3dw/fEQAEiNoCwIdIt2ia2Qwz225mO81sfpzzU8ys2sxKmn89GGUeAIinurZB1bUNvmMACAy1BYAPka3gmVmKpCclXSepTNIGM3vVObe13dTfOedujCoHAHRkwwcHJEnTRg3wnARASKgtAHyIcgVvkqSdzrlS59wxScsk3RThzwMAAACApBZlg5cjaU+r47LmsfYmm9l7ZvbvZjY6wjwAAAAAELQob7JiccZcu+N3JQ11zh02sy9IWiFpxAlvZDZL0ixJGjJkSBfHBAAAAIAwRLmCVyZpcKvjXEl7W09wzh1yzh1u/nqVpFQz69f+jZxzS5xzBc65gv79+0cYGQAAAAASV5QreBskjTCzPEnlkm6VNLP1BDO7WNI+55wzs0mKNZyVEWYCgBOMycnyHQFAgKgtAHyIrMFzzjWa2VxJqyWlSHrGObfFzGY3n18s6SuS/qeZNUqqlXSrc679Nk4AiNTFWWm+IwAIELUFgA+WaP1UQUGBKyoq8h0DQEAOHjkmSerds5vnJABCQm0BEBUzK3bOFcQ7F+mDzgEgERTvPqji3Qd9xwAQGGoLAB9o8AAAAAAgEDR4AAAAABAIGjwAAAAACAQNHgAAAAAEIsrn4AFAQhg/ONt3BAABorYA8IEGD0DS65/R3XcEAAGitgDwgS2aAJJeRU29KmrqfccAEBhqCwAfaPAAJL339lTpvT1VvmMACAy1BYAPNHgAAAAAEAgaPAAAAAAIBA0eAAAAAASCBg8AAAAAAsFjEgAkvfyhvX1HABAgagsAH2jwACS93j27+Y4AIEDUFgA+sEUTQNL7uLpOH1fX+Y4BIDDUFgA+sIIHIOn9sbxaknRxVprnJABCQm0B4AMreAAAAAAQCBo8AAAAAAgEDR4AAAAABIIGDwAAAAACwU1WACS9y/L6+I4AIEDUFgA+0OABSHpZ6am+IwAIELUFgA9s0QSQ9MoOHlXZwaO+YwAIDLUFgA+s4AFIeu9/VCNJyu3dw3MSACGhtgDwgRU8AAAAAAgEDR4AAAAABIIGDwAAAAACQYMHAAAAAIGItMEzsxlmtt3MdprZ/FPMu8zMmszsK1HmAYB4Jn+mryZ/pq/vGAACQ20B4ENkd9E0sxRJT0q6TlKZpA1m9qpzbmuceY9JWh1VFgA4lZ7duaEwgK5HbQHgQ5QreJMk7XTOlTrnjklaJummOPPukfSSpE8izAIAJ7W78oh2Vx7xHQNAYKgtAHyIssHLkbSn1XFZ81gLM8uR9CVJiyPMAQCntGPfYe3Yd9h3DACBobYA8CHKBs/ijLl2xz+T9H3nXNMp38hslpkVmVlRRUVFV+UDAAAAgKBEuTm8TNLgVse5kva2m1MgaZmZSVI/SV8ws0bn3IrWk5xzSyQtkaSCgoL2TSIAAAAAQNE2eBskjTCzPEnlkm6VNLP1BOdc3qdfm9mzkla2b+4AAAAAAJ0TWYPnnGs0s7mK3R0zRdIzzrktZja7+TzX3QEAAABAFzLnEmvHY0FBgSsqKvIdA0BA6hpilwGnpaZ4TgIgJNQWAFExs2LnXEG8czygBUDSS/YPX7c8vU6S9MLdkz0nAcKS7LUFgB9R3kUTABJCacVhlVZwK3MAXYvaAsAHGjwASa+04ohKK3gYMYCuRW0B4AMNHgAAAAAEggYPAAAAAAJBgwcAAAAAgaDBAwAAAIBA8JgEAElvysj+viMACBC1BYAPNHgAkt6FKWxmAND1qC0AfKDyAEh6O/bVaMe+Gt8xAASG2gLABxo8AElvd+VR7a486jsGgMBQWwD4QIMHAAAAAIGgwQMAAACAQNDgAQAAAEAgaPAAAAAAIBA8JgFA0ps2aoDvCAACRG0B4AMreAAAAAAQCBo8AElv20eHtO2jQ75jAAgMtQWADzR4AJJe+cFalR+s9R0DQGCoLQB8oMEDAAAAgEDQ4AEAAABAIGjwAAAAACAQPCYBQNJLSTHfEQAEiNoCwAcaPABJb+rIi3xHABAgagsAH9iiCQBJbMXGcm38sEp/+OCArnj0La3YWO47EgAAOAs0eACS3h/Lq/XH8mrfMc65FRvLtWD5Zh1rOi5JKq+q1YLlm2nygC6SrLUFgF+davDM7P92ZgwAEtHH1XX6uLrOd4xzbuHq7aptaGozVtvQpIWrt3tKBIQlWWsLAL86u4I3uvWBmaVIyu/6OACAc2VvVfwHMJ9sHAAAnP9O2eCZ2QIzq5E0zswONf+qkfSJpFfOSUIAQCQGZaef1jgAADj/nbLBc8792DmXIWmhcy6z+VeGc66vc25BR29uZjPMbLuZ7TSz+XHO32Rmm8ysxMyKzOzKs/hvAQCchnnTRyo9NaXNWHpqiuZNH+kpEQAAOFud3aK53syyPj0ws2wzKzzVNzRv43xS0g2SRkn6OzMb1W7am5LGO+cmSLpT0i86mQcAukz31AvUPTX57jlVODFHP755rLqlxP7bc7LT9eObx6pwYo7nZDif3fL0Ot3y9DrfMRJCstYWAH519jl4DznnXv70wDlXZWYPSVpxiu+ZJGmnc65UksxsmaSbJG1t9T6HW83vKcl1Mg8AdJmrRvT3HcGbwok5en79h5KkF+6e7DkNEJZkri0A/OnsPyvFm9dRc5gjaU+r47LmsTbM7Etm9r6k1xRbxQMAAAAAnIHOruAVmdnjim25dJLukVTcwfdYnLETVuiaVwZfNrOrJf1Q0rQT3shslqRZkjRkyJBORgaAzinZUyVJmjA422uO80pZsVSyVKraLWUPlSbMlHK5eTJwOqgtAHzo7ArePZKOSXpB0q8l1Uqa08H3lEka3Oo4V9Lek012zv2XpM+YWb8455Y45wqccwX9+7PdAUDX2l9Tr/019b5jnD/KiqU1j0i1lVLmwNjrmkdi4wA6jdoCwIdOreA5545Imm9mvdpdN3cqGySNMLM8SeWSbpU0s/UEM/srSX92zjkz+6ykbpIqO50eAND1SpZKaZlSWvO9tT59LVnKKh4AAOe5Tq3gmdnnzWyrmm+QYmbjzeyfTvU9zrlGSXMlrZa0TdKvnXNbzGy2mc1unvZlSX80sxLFtn/e4pzjRisA4FPVbql7Rtux7hmxcQAAcF7r7DV4P5U0XdKrkuSce6/5mrlTcs6tkrSq3djiVl8/JumxTqcFAEQve2hsW2Za1l/G6mti4wAA4LzW6YezOOf2tBtq6uIsAOBFj24p6tEtpeOJyWLCTKnukFRXLbnjsde6Q7FxAJ1GbQHgQ2dX8PaY2eclOTPrJulexbZdAkDC+/xfnXBvp+SWmy9Nva/tXTQn38P1d8BporYA8KGzDd5sSf+o2HPsyiT9pzq+iyYAIFHl5tPQAQCQgE7Z4JnZY86570ua6py77RxlAoBzqnj3AUlS/tA+npMACAm1BYAPHV2D9wUzS5W04FyEAQAfDh5p0MEjDb5jAAgMtQWADx1t0fwPSfsl9TSzQ5JMkvv01TmXGXE+AAAAAEAndbSCd79zLkvSa865TOdcRuvXcxEQAAAAANA5HTV465pfD0UdBAAAAABwdjraotnNzG6X9Hkzu7n9Sefc8mhiAcC5k5HW2RsKA0DnUVsA+NBR5Zkt6TZJ2ZL+tt05J4kGD0DC+9zwvr4jAAgQtQWAD6ds8Jxzb0t628yKnHP/fI4yAQAAAADOwCmvwTOz70mSc+6fzeyr7c79KMpgAHCu/KG0Un8orfQdA0BgqC0AfOjoJiu3tvq6/bPwZnRxFgDwoqauUTV1jb5jAAgMtQWADx1dg2cn+TreMQAASDZlxVLJUqlqt5Q9VJowU8rN950KAJJWRyt47iRfxzsGAADJpKxYWvOIVFspZQ6Mva55JDYOAPCioxW88WZ2SLHVuvTmr9V8nBZpMgDAOfHC3ZN9R0CiKlkqpWVKaVmx409fS5ayigcAnnR0F82UcxUEAHzp3TPVdwQgIazYWK6NH1bpWNNxXfHoW3opc7suzs1rO6l7Rmy7JqgtALzgCZwAkl7+0D6+IwDnvRUby7Vg+WYdazouSSqvqtVbR9N09YWfKHfgxX+ZWF8TuxYP1BYAXnR0DR4AAIAWrt6u2oamNmPLGq5S2UcfS3XVkjsee607FLvRCgDACxo8AEnvnZ379c7O/b5jAOe1vVW1J4xtcn+lH9fdLKX3lQ59FHudeh/X3zWjtgDwgS2aAJLe0WNNHU8Cktyg7HSVx2ny9meNlW78Xx4Snf+oLQB8YAUPAAB0aN70kUpPbXvvtfTUFM2bPtJTIgBAPKzgAQCADhVOzJEkfe/FTTrWdFw52emaN31kyzgA4PxAgwcAADqlcGKOnl//oSSenwgA5ysaPABJr19Gd98RAASI2gLABxo8AElvwuBs3xEABIjaAsAHbrICAAAAAIGgwQOQ9H63o0K/21HhOwaAwFBbAPjAFk0ASa++4bjvCAACRG0B4EOkK3hmNsPMtpvZTjObH+f8bWa2qfnXO2Y2Pso8AAAAABCyyBo8M0uR9KSkGySNkvR3Zjaq3bQPJF3jnBsn6YeSlkSVBwAAAABCF+UK3iRJO51zpc65Y5KWSbqp9QTn3DvOuYPNh7+XlBthHgAAAAAIWpTX4OVI2tPquEzS504x/79L+vcI8wBAXBdnpfmOACBA1BYAPkTZ4FmcMRd3otlUxRq8K09yfpakWZI0ZMiQrsoHAJKkMTlZviMACBC1BYAPUW7RLJM0uNVxrqS97SeZ2ThJv5B0k3OuMt4bOeeWOOcKnHMF/fv3jyQsAAAAACS6KBu8DZJGmFmemXWTdKukV1tPMLMhkpZL+rpz7k8RZgGAk1qz/ROt2f6J7xgAAkNtAeBDZFs0nXONZjZX0mpJKZKecc5tMbPZzecXS3pQUl9J/2RmktTonCuIKhMAxNPUFHf3OACcFWoLAB8ifdC5c26VpFXtxha3+vouSXdFmQEAAAAAkkWkDzoHAAAAAJw7NHgAAAAAEIhIt2gCQCLI6Z3uOwKAAFFbAPhAgwcg6V06MNN3BAABorYA8IEtmgAAAAAQCBo8AEnvja379MbWfb5jAAgMtQWADzR4AAAAABAIrsEDAACd9sLdk31HAACcAit4AAAAABAIGjwAAAAACARbNAEkvaF9e/iOACBA1BYAPtDgAUh6IwZk+I4AIEDUFgA+sEUTQNJrbDquxqbjvmMACAy1BYAPNHgAkt7a7RVau73CdwwAgaG2APCBBg8I1C1Pr9MtT6/zHQMAAADnEA0eAAAAAASCBg8AAAAAAkGDBwAAAACB4DEJAJLe8P49fUcAEKBkri0rNpZr4ert2ltVq0HZ6Zo3faQKJ+ZIZcVSyVKpareUPVSaMFPKzfcdFwgKDR6ApDe8fy/fEQAEKFlry4qN5VqwfLNqG5okSeVVtVqwfLOyD2zSlL0/l9IypcyBUm2ltOYRaep9NHlAF2KLJoCkV9fQpLrmDyIA0FWStbYsXL29pbn7VG1DkyrfeTbW3KVlSXZB7DUtM7aiB6DL0OABSHpv79ivt3fs9x0DQGCStbbsraqNO97n2EdS94y2g90zYts1AXQZGjwAAAB0mUHZ6XHHD3QbKNXXtB2sr4ldiwegy9DgAQAAoMvMmz5S6akpbcbSU1PU9/PflOoOSXXVkjsee607FLvRCoAuw01WAAAA0GUKJ+ZI0gl30ZwyMUcqu6jtXTQn38MNVoAuRoOXQG55ep0k6YW7J3tOAgAAcHKFE3NaGr02cvNp6ICI0eABSHojBiTnrcwBRIvaAsAHGjwASW9o3+R9GDGA6FBbAPjATVYAJL0j9Y06Ut/oOwaAwFBbAPgQaYNnZjPMbLuZ7TSz+XHO/7WZrTOzejP7bpRZAOBk1v25Uuv+XOk7BoDAUFsA+BDZFk0zS5H0pKTrJJVJ2mBmrzrntraadkDSvZIKo8oBAAAAAMkiyhW8SZJ2OudKnXPHJC2TdFPrCc65T5xzGyQ1RJgDAAAAAJJClA1ejqQ9rY7LmscAAAAAABGIssGzOGPujN7IbJaZFZlZUUVFxVnGAgAAAIAwRfmYhDJJg1sd50raeyZv5JxbImmJJBUUFJxRkwgAJ/PXAzN8RwAQIGoLAB+ibPA2SBphZnmSyiXdKmlmhD8PAM5Ibu8eviMACBC1BYAPkTV4zrlGM5srabWkFEnPOOe2mNns5vOLzexiSUWSMiUdN7P/LWmUc+5QVLkAoL3q2th9nrLSUz0nARASagsAH6JcwZNzbpWkVe3GFrf6+mPFtm4CgDcbPjggSZo2aoDnJABCQm0B4EOkDzoHAAAAAJw7NHgAAAAAEAgaPAAAAAAIBA0eEKAVG8u18cMq/eGDA7ri0be0YmO570gAAAA4ByK9yQqAc2/FxnItWL5Zx5qOS5LKq2q1YPlmSVLhxByf0c5bY3KyfEcAECBqCwAfWMEDArNw9XbVNjS1GattaNLC1ds9JTr/XZyVpouz0nzHABAYagsAH2jwgMDsrao9rXFIB48c08Ejx3zHABAYagsAH2jwgMAMyk4/rXFIxbsPqnj3Qd8xAASG2gLABxo8IDDzpo9UempKm7H01BQ9clm9tPI70r9+JfZaVuwpIQAAAKJCg5foyor50I42Cifm6Mc3j1W3lNj/3jnZ6XpqitOUvT+XaiulzIGx1zWP8PcFAAAgMDR4iaysOPYhnQ/taKdwYo4mDsnW5/L66P/Nv1ZTat+Q0jKltCzJLoi9pmVKJUt9RwUAAEAXosFLZCVL+dCOzqnaLXXPaDvWPSM2DgAAgGDwHLxEVrU7tnLXGh/aEU/20NgKb1qrZzLV18TGofGDs31HABAgagsAH1jBS2TZQ2Mf0lvjQzvimTBTqjsk1VVL7njste5QbBzqn9Fd/TO6+44BIDDUFgA+0OAliBUby7Xxwyr94YMDuuLRt7RiYzkf2tF5ufnS1Puk9L7SoY9ir1Pvi41DFTX1qqip9x0DQGCoLQB8YItmAlixsVwLlm/WsabjkqTyqlotWL5ZunmsCqfeF7vmrmp3bOVu8j18aEd8ufn83TiJ9/ZUSZKmjRrgNwiAoFBbAPhAg5cAFq7ertqGpjZjtQ1NWrh6uwrnX8uHdgAAAACS2KKZEPZW1Z7WOAAAAIDkRIOXAAZlp5/WOAAAAIDkRIOXAOZNH6n01JQ2Y+mpKZo3faSnRAAAAADOR1yDlwAKJ+ZIkr734iYdazqunOx0zZs+smUcwNnJH9rbdwQAAaK2APCBBi9BFE7M0fPrP5QkvXD3ZM9pgLD07tnNdwQAAaK2APCBLZoAkt7H1XX6uLrOdwwAgaG2APCBFTwASe+P5dWSpIuz0jwnARASagsAH1jBAwAAAIBA0OABAAAAQCBo8AAAAAAgEDR4AAAAABCISBs8M5thZtvNbKeZzY9z3sxsUfP5TWb22SjzAEA8l+X10WV5fXzHABAYagsAHyK7i6aZpUh6UtJ1ksokbTCzV51zW1tNu0HSiOZfn5P0VPMrAJwzWempviMACBC1BYAPUa7gTZK00zlX6pw7JmmZpJvazblJ0i9dzO8lZZvZwAgzAcAJyg4eVdnBo75jAAgMtQWAD1E2eDmS9rQ6LmseO905ABCp9z+q0fsf1fiOASAw1BYAPkTZ4FmcMXcGc2Rms8ysyMyKKioquiQcAAAAAIQmygavTNLgVse5kvaewRw555Y45wqccwX9+/fv8qAAAAAAEIIoG7wNkkaYWZ6ZdZN0q6RX2815VdI3mu+mebmkaufcRxFmAgAAAIBgRXYXTedco5nNlbRaUoqkZ5xzW8xsdvP5xZJWSfqCpJ2Sjkq6I6o8AAAAABC6yBo8SXLOrVKsiWs9trjV107SnCgzAMnqhbsn+46QMCZ/pq/vCAACRG0B4EOkDR4AJIKe3SmFALoetQWAD1FegwcACWF35RHtrjziOwaAwFBbAPjAPy0BSHo79h2WJA3t29NzEgAhobYA8IEVPAAAAAAIBA0eAAAAAASCBg8AAAAAAkGDBwAAAACB4CYrAJLelSP6+Y4AIEDUFgA+0OABSHppqSm+IwAIELUFgA9s0QSQ9EorDqu04rDvGAACQ20B4AMNHoCkV1pxRKUVPIwYQNeitgDwgS2aCeSFuyf7jgAAAADgPMYKHgAAAAAEggYPAAAAAAJBgwcAAAAAgeAaPABJb8rI/r4jAAgQtQWADzR4AJLehSlsZgDQ9agtAHyg8gBIejv21WjHvhrfMQAEhtoCwAcaPABJb3flUe2uPOo7BoDAUFsA+ECDBwAAAACBoMEDAAAAgEDQ4AEAAABAIGjwAAAAACAQ5pzzneG0mFmFpN2+cyBh9JO033cIAMGhtgCIArUFnTXUORf3YZsJ1+ABp8PMipxzBb5zAAgLtQVAFKgt6Aps0QQAAACAQNDgAQAAAEAgaPAQuiW+AwAIErUFQBSoLThrXIMHAAAAAIFgBQ8AAAAAAkGDBwAAAACBoMEDAAAAgEDQ4CEoZpZiZv9oZlvMbLOZDfedCUDio7YAiAK1BVGgwUNoFkgqdc6NlrRI0rc85wEQBmoLgChQW9DlLvQdAOgqZtZT0pecc/nNQx9I+qLHSAACQG0BEAVqC6JCg4eQTJM02MxKmo/7SHrDXxwAgaC2AIgCtQWRYIsmQjJB0oPOuQnOuQmS/lNSiZn1NLPnzOznZnab14QAEtEExa8tw83sn83sRa/pACSqCYpfWwqbP7O8YmbXe02IhESDh5D0lnRUkszsQknXS/qNpJslveic+x+S/pu/eAASVNza4pwrdc79d6/JACSyk9WWFc2fWb4p6RZ/8ZCoaPAQkj9Jurz5629Les0594GkXEl7msebfAQDkNBOVlsA4Gx0VFvul/TkOU+FhEeDh5A8L+mzZrZT0jhJ/6d5vEyxJk/i7zyA03ey2gIAZyNubbGYxyT9u3PuXZ8BkZjMOec7AxCp5rtUPSGpTtLbzrlfeY4EIABm1lfSI5Kuk/QL59yPPUcCEAAzu1fS7ZI2SCpxzi32HAkJhgYPAAAAAALBdjUAAAAACAQNHgAAAAAEggYPAAAAAAJBgwcAAAAAgaDBAwAAAIBA0OABAAAAQCBo8AAAAAAgEDR4AAAAABAIGjwAAAAACMT/B1HSZiIhWdHSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Some plotting boilerplate code\n",
    "plt.figure(figsize=(15, 5))\n",
    "plt.errorbar(np.arange(n_periods*n_treatments)-.04, est.intercept_, yerr=(conf_ints[1] - est.intercept_,\n",
    "                                                    est.intercept_ - conf_ints[0]), fmt='o', label='DynamicDML')\n",
    "plt.errorbar(np.arange(n_periods*n_treatments), true_effect.flatten(), fmt='o', alpha=.6, label='Ground truth')\n",
    "for t in np.arange(1, n_periods):\n",
    "    plt.axvline(x=t * n_treatments - .5, linestyle='--', alpha=.4)\n",
    "plt.xticks([t * n_treatments - .5 + n_treatments/2 for t in range(n_periods)],\n",
    "           [\"$\\\\theta_{}$\".format(t) for t in range(n_periods)])\n",
    "plt.gca().set_xlim([-.5, n_periods*n_treatments - .5])\n",
    "plt.ylabel(\"Effect\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Example Usage with Heterogeneous Treatment Effects on Time-Invariant Unit Characteristics\n",
    "\n",
    "We can also estimate treatment effect heterogeneity with respect to the value of some subset of features $X$ in the initial period. Heterogeneity is currently only supported with respect to such initial state features. This for instance can support heterogeneity with respect to time-invariant unit characteristics. In that case you can simply pass as $X$ a repetition of some unit features that stay constant in all periods. You can also pass time-varying features, and their time varying component will be used as a time-varying control. However, heterogeneity will only be estimated with respect to the initial state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 DGP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define additional DGP parameters\n",
    "het_strength = .5\n",
    "het_inds = np.arange(n_x - n_treatments, n_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate data\n",
    "dgp = DynamicPanelDGP(n_periods, n_treatments, n_x).create_instance(\n",
    "            s_x, hetero_strength=het_strength, hetero_inds=het_inds, random_seed=12)\n",
    "Y, T, X, W, groups = dgp.observational_data(n_panels, s_t=s_t, random_seed=1)\n",
    "ate_effect = dgp.true_effect\n",
    "het_effect = dgp.true_hetero_effect[:, het_inds + 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Train Estimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "est = DynamicDML(\n",
    "    model_y=LassoCV(cv=3),\n",
    "    model_t=MultiTaskLassoCV(cv=3),\n",
    "    cv=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<econml.panel.dml._dml.DynamicDML at 0x19d2ae7e5e0>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est.fit(Y, T, X=X, W=W, groups=groups, inference=\"auto\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Coefficient Results</caption>\n",
       "<tr>\n",
       "       <td></td>       <th>point_estimate</th> <th>stderr</th>  <th>zstat</th>  <th>pvalue</th> <th>ci_lower</th> <th>ci_upper</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X0|(T0)$_0$</th>      <td>0.009</td>      <td>0.045</td>  <td>0.203</td>   <td>0.839</td>  <td>-0.065</td>    <td>0.083</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X0|(T1)$_0$</th>      <td>0.017</td>      <td>0.042</td>  <td>0.416</td>   <td>0.677</td>  <td>-0.051</td>    <td>0.086</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X0|(T0)$_1$</th>     <td>-0.001</td>      <td>0.041</td> <td>-0.035</td>   <td>0.972</td>  <td>-0.069</td>    <td>0.067</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X0|(T1)$_1$</th>     <td>-0.031</td>      <td>0.041</td>  <td>-0.76</td>   <td>0.447</td>  <td>-0.099</td>    <td>0.036</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X0|(T0)$_2$</th>     <td>-0.306</td>      <td>0.008</td> <td>-36.667</td>   <td>0.0</td>    <td>-0.32</td>   <td>-0.292</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X0|(T1)$_2$</th>      <td>0.158</td>      <td>0.008</td> <td>19.656</td>    <td>0.0</td>    <td>0.145</td>    <td>0.171</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X1|(T0)$_0$</th>      <td>0.017</td>      <td>0.044</td>  <td>0.378</td>   <td>0.706</td>  <td>-0.056</td>    <td>0.09</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X1|(T1)$_0$</th>     <td>-0.007</td>      <td>0.045</td> <td>-0.164</td>   <td>0.87</td>   <td>-0.082</td>    <td>0.067</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X1|(T0)$_1$</th>     <td>-0.034</td>      <td>0.042</td> <td>-0.821</td>   <td>0.412</td>  <td>-0.103</td>    <td>0.034</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X1|(T1)$_1$</th>     <td>-0.025</td>      <td>0.042</td>  <td>-0.6</td>    <td>0.549</td>  <td>-0.095</td>    <td>0.044</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X1|(T0)$_2$</th>     <td>-0.302</td>      <td>0.008</td> <td>-35.72</td>    <td>0.0</td>   <td>-0.316</td>   <td>-0.288</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>X1|(T1)$_2$</th>      <td>0.156</td>      <td>0.008</td> <td>18.801</td>    <td>0.0</td>    <td>0.142</td>    <td>0.169</td> \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<caption>CATE Intercept Results</caption>\n",
       "<tr>\n",
       "             <td></td>             <th>point_estimate</th> <th>stderr</th>   <th>zstat</th>  <th>pvalue</th> <th>ci_lower</th> <th>ci_upper</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T0)$_0$</th>      <td>0.024</td>      <td>0.036</td>   <td>0.653</td>   <td>0.513</td>  <td>-0.036</td>    <td>0.084</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T1)$_0$</th>     <td>-0.033</td>      <td>0.036</td>  <td>-0.929</td>   <td>0.353</td>  <td>-0.092</td>    <td>0.025</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T0)$_1$</th>     <td>-0.105</td>      <td>0.034</td>  <td>-3.067</td>   <td>0.002</td>  <td>-0.162</td>   <td>-0.049</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T1)$_1$</th>      <td>0.037</td>      <td>0.034</td>   <td>1.079</td>   <td>0.281</td>  <td>-0.019</td>    <td>0.093</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T0)$_2$</th>     <td>-0.743</td>      <td>0.005</td> <td>-140.503</td>   <td>0.0</td>   <td>-0.752</td>   <td>-0.734</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cate_intercept|(T1)$_2$</th>      <td>0.48</td>       <td>0.005</td>  <td>91.061</td>    <td>0.0</td>    <td>0.472</td>    <td>0.489</td> \n",
       "</tr>\n",
       "</table><br/><br/><sub>A linear parametric conditional average treatment effect (CATE) model was fitted:<br/>$Y = \\Theta(X)\\cdot T + g(X, W) + \\epsilon$<br/>where for every outcome $i$ and treatment $j$ the CATE $\\Theta_{ij}(X)$ has the form:<br/>$\\Theta_{ij}(X) = \\phi(X)' coef_{ij} + cate\\_intercept_{ij}$<br/>where $\\phi(X)$ is the output of the `featurizer` or $X$ if `featurizer`=None. Coefficient Results table portrays the $coef_{ij}$ parameter vector for each outcome $i$ and treatment $j$. Intercept Results table portrays the $cate\\_intercept_{ij}$ parameter.</sub>"
      ],
      "text/plain": [
       "<class 'econml.utilities.Summary'>\n",
       "\"\"\"\n",
       "                       Coefficient Results                        \n",
       "==================================================================\n",
       "            point_estimate stderr  zstat  pvalue ci_lower ci_upper\n",
       "------------------------------------------------------------------\n",
       "X0|(T0)$_0$          0.009  0.045   0.203  0.839   -0.065    0.083\n",
       "X0|(T1)$_0$          0.017  0.042   0.416  0.677   -0.051    0.086\n",
       "X0|(T0)$_1$         -0.001  0.041  -0.035  0.972   -0.069    0.067\n",
       "X0|(T1)$_1$         -0.031  0.041   -0.76  0.447   -0.099    0.036\n",
       "X0|(T0)$_2$         -0.306  0.008 -36.667    0.0    -0.32   -0.292\n",
       "X0|(T1)$_2$          0.158  0.008  19.656    0.0    0.145    0.171\n",
       "X1|(T0)$_0$          0.017  0.044   0.378  0.706   -0.056     0.09\n",
       "X1|(T1)$_0$         -0.007  0.045  -0.164   0.87   -0.082    0.067\n",
       "X1|(T0)$_1$         -0.034  0.042  -0.821  0.412   -0.103    0.034\n",
       "X1|(T1)$_1$         -0.025  0.042    -0.6  0.549   -0.095    0.044\n",
       "X1|(T0)$_2$         -0.302  0.008  -35.72    0.0   -0.316   -0.288\n",
       "X1|(T1)$_2$          0.156  0.008  18.801    0.0    0.142    0.169\n",
       "                             CATE Intercept Results                            \n",
       "===============================================================================\n",
       "                        point_estimate stderr  zstat   pvalue ci_lower ci_upper\n",
       "-------------------------------------------------------------------------------\n",
       "cate_intercept|(T0)$_0$          0.024  0.036    0.653  0.513   -0.036    0.084\n",
       "cate_intercept|(T1)$_0$         -0.033  0.036   -0.929  0.353   -0.092    0.025\n",
       "cate_intercept|(T0)$_1$         -0.105  0.034   -3.067  0.002   -0.162   -0.049\n",
       "cate_intercept|(T1)$_1$          0.037  0.034    1.079  0.281   -0.019    0.093\n",
       "cate_intercept|(T0)$_2$         -0.743  0.005 -140.503    0.0   -0.752   -0.734\n",
       "cate_intercept|(T1)$_2$           0.48  0.005   91.061    0.0    0.472    0.489\n",
       "-------------------------------------------------------------------------------\n",
       "\n",
       "<sub>A linear parametric conditional average treatment effect (CATE) model was fitted:\n",
       "$Y = \\Theta(X)\\cdot T + g(X, W) + \\epsilon$\n",
       "where for every outcome $i$ and treatment $j$ the CATE $\\Theta_{ij}(X)$ has the form:\n",
       "$\\Theta_{ij}(X) = \\phi(X)' coef_{ij} + cate\\_intercept_{ij}$\n",
       "where $\\phi(X)$ is the output of the `featurizer` or $X$ if `featurizer`=None. Coefficient Results table portrays the $coef_{ij}$ parameter vector for each outcome $i$ and treatment $j$. Intercept Results table portrays the $cate\\_intercept_{ij}$ parameter.</sub>\n",
       "\"\"\""
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average effect of default policy:-0.42\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "A scalar was specified but there are multiple treatments; the same value will be used for each treatment.  Consider specifyingall treatments, or using the const_marginal_effect method.\n",
      "A scalar was specified but there are multiple treatments; the same value will be used for each treatment.  Consider specifyingall treatments, or using the const_marginal_effect method.\n"
     ]
    }
   ],
   "source": [
    "# Average treatment effect for test points\n",
    "X_test = X[np.arange(0, 25, 3)]\n",
    "print(f\"Average effect of default policy:{est.ate(X=X_test):0.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Effect of target policy over baseline policy for test set:\n",
      " [-0.37368525 -0.30896804 -0.43030363 -0.52252401 -0.42849622 -0.48790877\n",
      " -0.34417987 -0.51804937 -0.36806744]\n"
     ]
    }
   ],
   "source": [
    "# Effect of target policy over baseline policy\n",
    "# Must specify a treatment for each period\n",
    "baseline_policy = np.zeros((1, n_periods * n_treatments))\n",
    "target_policy = np.ones((1, n_periods * n_treatments))\n",
    "eff = est.effect(X=X_test, T0=baseline_policy, T1=target_policy)\n",
    "print(\"Effect of target policy over baseline policy for test set:\\n\", eff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0.02374269, -0.03302781, -0.10526464,  0.03675719, -0.74294675,\n",
       "         0.48025068]),\n",
       " array([[ 0.00914226,  0.01675409],\n",
       "        [ 0.01732804, -0.00741467],\n",
       "        [-0.00143705, -0.03431712],\n",
       "        [-0.03136295, -0.02536834],\n",
       "        [-0.30581311, -0.30189654],\n",
       "        [ 0.15773252,  0.15564665]]))"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Coefficients: intercept is of shape n_treatments*n_periods\n",
    "# coef_ is of shape (n_treatments*n_periods, n_hetero_inds).\n",
    "# first n_treatment rows are from first period, next n_treatment\n",
    "# from second period, etc.\n",
    "est.intercept_, est.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Confidence intervals\n",
    "conf_ints_intercept = est.intercept__interval(alpha=0.05)\n",
    "conf_ints_coef = est.coef__interval(alpha=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Performance Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parse true parameters in array of shape (n_treatments*n_periods, 1 + n_hetero_inds)\n",
    "# first column is the intercept\n",
    "true_effect_inds = []\n",
    "for t in range(n_treatments):\n",
    "    true_effect_inds += [t * (1 + n_x)] + (list(t * (1 + n_x) + 1 + het_inds) if len(het_inds)>0 else [])\n",
    "true_effect_params = dgp.true_hetero_effect[:, true_effect_inds]\n",
    "true_effect_params = true_effect_params.reshape((n_treatments*n_periods, 1 + het_inds.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# concatenating intercept and coef_\n",
    "param_hat = np.hstack([est.intercept_.reshape(-1, 1), est.coef_])\n",
    "lower = np.hstack([conf_ints_intercept[0].reshape(-1, 1), conf_ints_coef[0]])\n",
    "upper = np.hstack([conf_ints_intercept[1].reshape(-1, 1), conf_ints_coef[1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAFgCAYAAABNIolGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAySklEQVR4nO3de7hdZX0v+u+bCwmXJIRLkBAwZFc93EIskepuqyAXQauwt9tocXvBosdHKPt42CK7IiBlt9pa3c2RVqkX8HjBFKqo6HZD1BZOqRhsGmSjYDFCACEEiCAkJOQ9f8xFClnJzGWurDHWnJ/P84xnzjHmWHP8Ruab75rzt8YYs9RaAwAAANBm45ouAAAAAGBrNDAAAACA1tPAAAAAAFpPAwMAAABoPQ0MAAAAoPUmNF1AN/vss0+dPXt202XAyFu3rnM7cWKzdUCShx5fmyTZZ49JDVcCQ2QkLSMnaR05SZ+75ZZbHqq17rvp8lY3MGbPnp0lS5Y0XQZAX/ubf7grSfLOl89puBKAdpKTAKOrlPKLzS13CgkAAADQehoY0IR//dfOBMBwMhKgOznJgGr1KSTQt55+uukKANpLRgJ0JycZUBoYAAAA9L1169ZlxYoVWbNmTdOlMGTy5MmZNWtWJm7jBWk1MAAAAOh7K1asyJQpUzJ79uyUUpouZ+DVWrNq1aqsWLEiBx988Db9jGtgAAAA0PfWrFmTvffeW/OiJUop2XvvvbfriBhHYEATpkxpugKA9pKRAN3JyR2medEu2/t6aGBAE/bfv+kKANpLRgJ0JycZUE4hAQAAgFHwwAMP5LTTTsucOXNy1FFH5WUve1m++tWvjmoNy5cvz+GHH/6cZbfeemvmzZuXefPmZa+99srBBx+cefPm5fjjj9/m5/zSl760cf7yyy/PWWedNaJ1JxoY0Iw77+xMAAwnIwG6k5Oj5sFfrcmCT92UBx/r/ZtLaq059dRT8/KXvzx33XVXbrnlllx55ZVZsWLFsHXXr1/f8/a2xxFHHJGlS5dm6dKled3rXpc///M/z9KlS3P99ddvU02bNjB2Fg0MaEKtnQmA4WQkQHejlZPLFiUfPzy5aM/O7bJFO3+bLbNw8Z354fKHs/D63htG3/3ud7PLLrvk3e9+98Zlz3/+8/OHf/iHSTpHLbzhDW/Ia1/72px44ol5+OGHc+qpp2bu3Ll56UtfmmXLliVJLrroonz0ox/d+ByHH354li9fnuXLl+eQQw7JO9/5zhx22GE58cQT8+STTyZJbrnllhx55JF52ctelksvvXSbaz7mmGPyR3/0R3nFK16Rv/zLv8zb3/72XHXVVRsf32OPPZIk5513Xm644YbMmzcvH//4x5Mk9913X0466aS84AUvyLnnnruD/2rPpYEBAADAcy1blHzj7GT1PUlq5/YbZw9ME+NF5387s8+7Nl/4wd2pNfnCD+7O7POuzYvO//YOP+dtt92W3/zN3+y6zk033ZQrrrgi3/3ud3PhhRfmxS9+cZYtW5Y/+ZM/yVvf+tatbuPOO+/MmWeemdtuuy177rlnrr766iTJ6aefnoULF+amm27a7rofffTR/P3f/33OOeecLa7z4Q9/OL/7u7+bpUuX5r3vfW+SZOnSpfnKV76SW2+9NV/5yldyzz33bPe2N6WBAQAAwHMtvjhZ9+Rzl617srN8ANxw7rF53byZmTyx85F58sRxOWXezNzw/mNHbBtnnnlmjjzyyLzkJS/ZuOyEE07IXnvtlSS58cYb85a3vCVJ8spXvjKrVq3K6tWruz7nM9euSJKjjjoqy5cvz+rVq/Poo4/mFa94RZJsfM5t9cY3vnG71n/Gcccdl2nTpmXy5Mk59NBD84tf/GKHnufZNDAAAAB4rtXDr8vQdXmfmTF1cqZMmpC16zdk0oRxWbt+Q6ZMmpAZUybv8HMedthh+dGPfrRx/tJLL83ixYuzcuXKjct23333jffrZk4TKqVkwoQJ2bBhw8Zla9b82/U5Jk2atPH++PHjs379+tRae/r62GfX9Oxt11rz1FNPbfHnNldLrzQwoAnTpnUmAIaTkQDdjUJO3lv33q7l/eihx9fmzb/1/Hz1Pb+dN//W87Py8bU9Pd8rX/nKrFmzJn/913+9cdkTTzyxxfVf/vKX54tf/GKS5Pvf/3722WefTJ06NbNnz97YCPnRj36Un//85123u+eee2batGm58cYbk2Tjc+6I2bNn55ZbbkmSXHPNNVm3bl2SZMqUKXnsscd2+Hm31YSdvgVguP32a7oCgPaSkQDdjUJO7vHqi/PU/3xvdqn/9qH9qTIpe7x6ME4hSZJPvWX+xvuXnHp4lzW3TSklX/va1/Le9743f/Znf5Z99903u+++ez7ykY9sdv2LLroop59+eubOnZvddtstV1xxRZLk9a9/fT7/+c9n3rx5eclLXpIXvvCFW9325z73ubzjHe/Ibrvtlle96lU7vA/vfOc7c8opp+Too4/Occcdt/HojLlz52bChAk58sgj8/a3vz3Tp0/f4W10UzZ3WEpbzJ8/vy5ZsqTpMgD62t/8w11Jkne+fE7DlQC0k5xkUC36zF/k39/9V5mZVbkve+cfD3pPFvzBli/k2Ha33357DjnkkKbLYBObe11KKbfUWudvuq4jMKAJd9zRud2GbinAwJGRAN2NUk4u3uUVWfbiV+W0ow/Kl26+OysfW5MFO3WL0J0GBgAAAMOM9CkU0CsX8QQAAABaTwMDAAAAaD0NDAAAAKD1XAMDmrCTvlYIoC/ISIDu5OTo+dxrOrenX9tsHSRxBAY0Y999OxMAw8lIgO7k5Jg1fvz4zJs3L4cddliOPPLIfOxjH8uGDRsaqWXJkiU5++yzu64ze/bsHHHEETniiCNy6KGH5vzzz8/atWuTJMuXL08pJR/84Ac3rv/QQw9l4sSJOeuss5IkF110UT760Y+OWM0aGNCEDRs6EwDDyUiA7uTk6Fi2KFnxw+QXNyYfP7wz36Ndd901S5cuzW233Zbrrrsu3/rWt/KhD31oBIrdfvPnz8/ChQu3ut73vve93Hrrrbn55ptz11135V3vetfGx+bMmZNvfvObG+f/9m//NocddthOqTcZoQZGKeWkUspPSyk/K6Wc12W9l5RSni6l/KeR2C6MWT/7WWcCYDgZCdCdnNz5li1KvnF28nTnaIOsvqczPwJNjGfMmDEjl112WT7xiU+k1prf/d3fzdKlSzc+/tu//dtZtmxZLrroorzjHe/IMccckzlz5jyn6XDqqafmqKOOymGHHZbLLrts4/I99tgj73//+3PUUUfl+OOPz80337zx57/+9a8nSb7//e/n937v95Ikjz/+eE4//fQcccQRmTt3bq6++uph9e6xxx755Cc/ma997Wt5+OGHk3QaMoccckiWLFmSJPnKV76SBQsWjNi/0aZ6bmCUUsYnuTTJyUkOTfL7pZRDt7DeR5J8p9dtAgAAwE6z+OJk3ZPPXbbuyc7yETRnzpxs2LAhDz74YM4444xcfvnlSZI77rgja9euzdy5c5MkP/nJT/Kd73wnN998cz70oQ9l3bp1SZLPfvazueWWW7JkyZIsXLgwq1atSpL8+te/zjHHHJNbbrklU6ZMyfnnn5/rrrsuX/3qV3PBBRcMq+OP//iPM23atNx6661ZtmxZXvnKV2623qlTp+bggw/OnXfeuXHZm970plx55ZVZsWJFxo8fn5kzZ47kP9FzjMQRGEcn+Vmt9a5a61NJrkxyymbW+8MkVyd5cAS2CQAAADvH6hXbt7wHtdYkyRve8IZ885vfzLp16/LZz342b3/72zeu85rXvCaTJk3KPvvskxkzZuSBBx5IkixcuDBHHnlkXvrSl+aee+7Z2FjYZZddctJJJyVJjjjiiLziFa/IxIkTc8QRR2T58uXDarj++utz5plnbpyf3uVCsc/U+4yTTjop1113Xb785S/njW984w79G2yrkWhgHJDknmfNrxhatlEp5YAk/yHJJ0dgewAAALDzTJu1fct30F133ZXx48dnxowZ2W233XLCCSfkmmuuyaJFi3LaaadtXG/SpEkb748fPz7r16/P97///Vx//fW56aab8i//8i958YtfnDVr1iRJJk6cmFJKkmTcuHEbf37cuHFZv379sDpqrRvX7+axxx7L8uXL88IXvnDjsl122SVHHXVU/uIv/iKvf/3rd+wfYhuNRANjc3tZN5n/H0neX2t9eqtPVsq7SilLSilLVq5cOQLlAQAAwHY47oJk4q7PXTZx187yEbJy5cq8+93vzllnnbWxeXDGGWfk7LPPzkte8pLstddeXX9+9erVmT59enbbbbf85Cc/yT/90z/tcC0nnnhiPvGJT2ycf+SRR4at8/jjj+c973lPTj311GFHaJxzzjn5yEc+kr333nuHa9gWI9HAWJHkwGfNz0py3ybrzE9yZSlleZL/lOSvSimnbu7Jaq2X1Vrn11rn7+urgehXe+/dmQAYTkYCdCcnd765C5LXLkzGDx35MO3Azvzc3i5Q+eSTT278GtXjjz8+J554Yi688MKNjx911FGZOnVqTj/99K0+10knnZT169dn7ty5+eAHP5iXvvSlO1zX+eefn0ceeSSHH354jjzyyHzve9/b+Nixxx6bww8/PEcffXQOOuigfOpTnxr284cddlje9ra3bfa5L7nkksyaNWvj1Iuy6fkr2/0EpUxIckeS45Lcm+SHSU6rtd62hfUvT/LNWutVW3vu+fPn12euZgrAzvE3/3BXkuSdL5/TcCUA7SQnoT/cfvvtOeSQQ7bvhz73ms7t6deOfEGbcd999+WYY47JT37yk4wbNyJfGtp6m3tdSim31Frnb7ruhF43VmtdX0o5K51vFxmf5LO11ttKKe8eetx1L2BTz5x3NqHn/4IA/UdGAnQnJ0fPKDUukuTzn/98PvCBD+RjH/vYwDQvtteIjPha67eSfGuTZZttXNRa3z4S24Qx7a7OX3LyrIvfADBERgJ0Jyf70lvf+ta89a1vbbqMVtPWAQAAYCD0egkFRtb2vh4aGAAAAPS9yZMnZ9WqVZoYLVFrzapVqzJ58uRt/hknTQEAAND3Zs2alRUrVmTlypVNl8KQyZMnb9c3k2hgAAAA0PcmTpyYgw8+uOky6IEGBjRh332brgCgvWQkQHdykgGlgQFNmD696QoA2ktGAnQnJxlQLuIJTXjqqc4EwHAyEqA7OcmA0sCAJixf3pkAGE5GAnQnJxlQGhgAAABA62lgAAAAAK2ngQEAAAC0ngYGAAAA0Hq+RhWasN9+TVcA0F4yEqA7OcmA0sCAJkyb1nQFAO0lIwG6k5MMKKeQQBPWrOlMAAwnIwG6k5MMKA0MaMLdd3cmAIaTkQDdyUkGlAYGAAAA0HoaGAAAAEDraWAAAAAAraeBAQAAALSer1GFJjzveU1XANBeMhKgOznJgNLAgCZMndp0BQDtJSMBupOTDCinkEATnniiMwEwnIwE6E5OMqA0MKAJK1Z0JgCGk5EA3clJBpQGBgAAANB6GhgAAABA62lgAAAAAK2ngQEAAAC0nq9RhSbMnNl0BQDtJSMBupOTDCgNDGjCHns0XQFAe8lIgO7kJAPKKSTQhMcf70wADCcjAbqTkwyoEWlglFJOKqX8tJTys1LKeZt5/M2llGVD0z+WUo4cie3CmHXffZ0JgOFkJEB3cpIB1XMDo5QyPsmlSU5OcmiS3y+lHLrJaj9P8opa69wkf5zksl63CwAAAAyOkTgC4+gkP6u13lVrfSrJlUlOefYKtdZ/rLU+MjT7T0lmjcB2AQAAgAExEg2MA5Lc86z5FUPLtuQPknx7Sw+WUt5VSllSSlmycuXKESgPAAAAGOtGooFRNrOsbnbFUo5Np4Hx/i09Wa31slrr/Frr/H333XcEygMAAADGupH4GtUVSQ581vysJMOuKFNKmZvk00lOrrWuGoHtwtg1y1lUAFskIwG6k5MMqJE4AuOHSV5QSjm4lLJLkjcl+fqzVyilHJTk75K8pdZ6xwhsE8a23XbrTAAMJyMBupOTI+LBX63Jgk/dlAcfW9N0KSOuX/et5wZGrXV9krOSfCfJ7UkW1VpvK6W8u5Ty7qHVLkiyd5K/KqUsLaUs6XW7MKb96ledCYDhZCRAd3JyRCxcfGd+uPzhLLz+zqZLGXH9um8jcQpJaq3fSvKtTZZ98ln3z0hyxkhsC/rCL3/ZuZ06tdk6ANpIRgJ0Jyd78qLzv51XbfiHnDthUS7e5aHc98/75OwfLsh3xr08P73k5KbL68mLzv921q7fsHH+Cz+4O1/4wd2ZNGHcmN+3ZIQaGAAAADAW3Pzah7Pr//xMdqlrkySzykP56KTP5I9PPrzhynp3w7nH5ttfXpjj7/tk9s+q3J99snjm/5mTTju76dJGxEhcAwMAAADGhGn/+KcbmxfP2KWuzbT/708bqmjkzFj+9fz+Lz+aA8qqjCvJAeWhvOmXH82Mn3996z88BmhgAAAAMDhWr9i+5WPJ4os325zJ4osbKmhkaWAAAAAwOKZt4Wtot7R8LOnn5kw0MKAZBx3UmQAYTkYCdCcne3PcBcnEXZ+7bOKuneVjXT83Z6KBAc2YPLkzATCcjAToTk72Zu6C5LULk/GTOvPTDuzMz13QbF0joZ+bM/EtJNCM1as7t9OmNVsHQBvJSIDu5GTv5i7oj4bFpp7Zp8UXd04bmTar07zok33VwIAmPPBA59YvHYDhZCRAd3KSbvq1OROnkAAAAABjgAYGAAAA0HoaGAAAAEDraWAAAAAArecintCE2bObrgCgvWQkQHdykgGlgQFJ8rnXdG5Pv3Z0trfLLqOzHYCxaLQzcrR/BwD0yntJBpRTSNh2n3vNv73JozePPNKZ6I0xCf1JRgJ0JycZUI7AGGn+isO2WLmyczt9erN1ALSRjAToTk4yoByBAQAAALSeBgZAGzk9BgAAnkMDAwBGisYTAMBOo4EBAAAAtJ6LeEIT5sxpugKA9pKRAN3JSQaUBgY0YYL/egBbJCMBupOTDCinkEATVq3qTAAMJyPZGtebYdDJSQaUBgY0wS8dgC0bzYxctihZ8cPkFzcmHz+8Mw/Qdt5LMqA0MACArevHv3gvW5R84+zk6bWd+dX3dOb7pYnRj68ZAANNAwMAGEyLL07WPfncZeue7CwHAFpHA2MkOQyVNvIXOIDNW71i+5YDAI3SwBgp/X4Yaj/TeBqbvG5Ar6bN2r7lAECjNDBGisNQx6amGk+/8RudiR2jYQj9bbQy8rgLkom7PnfZxF07ywHazHtJBpQGxkjp98NQ+/Wv3U01nsaN60zsGA1D6G+jlZFzFySvXZinx0/OhiRPT5mVvHZhZzlAm3kvyYAakVFfSjmplPLTUsrPSinnbebxUkpZOPT4slLKb47Edlulnw9D7ee/djfVeFq5sjOxY/q9YQiDbjQzcu6CXDjhvfl3a76YC+dc2VfNiwfX75YFd5+SBx9b03QpI66f9w22ifeSDKieGxillPFJLk1ycpJDk/x+KeXQTVY7OckLhqZ3JfnrXrfbOv18GGo//7W7qcbTI490JnZMPzcMh3hzPgb165FqTRiljHzR+d/O7POuzRdWH5Gaki/84O7MPu/avOj8b+/0bY+GhQ8dlR8+OTMLr7+z6VJGXD/vG2wT7yUZUCNxBMbRSX5Wa72r1vpUkiuTnLLJOqck+Xzt+Kcke5ZS9h+BbbfGixZNydm/Pj0rNuyTDbVkxYZ9cvavT8+LFk1purSebXh083/V3tLyseSch0/JE3WX5yx7ou6Scx7edAjTKv3cMEySZYsy8d5/ypUb3peJC+f6IDwW9PORan3shnOPzevmzczksi5JMnniuJwyb2ZueP+xDVfWm35uzPTzvgGwdSPRwDggyT3Pml8xtGx71xnTbjj32GTughy/7i8yZ+0Xc3z9RMrcBWP+TVCS1Kmbf6m2tHwsef/7Ppi/nfm+3Fv3zoZacm/dJ1fNfF/ef+4Hmy5txPTjX/L7uWF4zgc/kCeuPjPT86uMK8n0dQ/kiavPzDkf/EDTpdFNPx+p1sdmTJ2cKZMmZG2dkEllfdau35ApkyZkxpTJTZfWk35tzCT9vW8AbN1INDDKZpbVHVins2Ip7yqlLCmlLFk5hs7r6tc3QUky/oQL81SZ9JxlT5VJGX/ChQ1VNHJmTJ2cO2acnN9ZuzCHPHVFfuephbljxsl98bo9ox8Ps+3nhuGf7fm17Faees6y3cpT+bM9v9ZMQWwb12UZsx56fG3ePO3H+epBV+fNv/X8rHx8bdMl9ayf35P0874BsHUTRuA5ViQ58Fnzs5LctwPrJElqrZcluSxJ5s+fv9kmR1s98ybotD3/d760/7lZ2S9/8Z67IJ+/8a68/sFLs2d+lUcn7perp78jZ/TJhc769XV70fnfztr1G5IckST5wg/uzhd+cHcmTRiXn15ycrPF9Wj4G9gJffMGdvxj927Xclpi2qzOaSObW06rfeot85PPdRryl5x6eMPVjJx+/d2W9Pe+AdDdSDQwfpjkBaWUg5Pcm+RNSU7bZJ2vJzmrlHJlkt9KsrrWev8IbLtV+vVNUJKc8Z7zks/dkCSZfvq1OaPhekZSI6/bC1+40zdxw7nH5pJv3Z7/9S+/yJo6MZMnjsurDntePvCaQ3b6tkdD376B9UF4bDrugs41L559Gkk/XZdltI1CRva7fn5P8qkj70qWfyx5YG0uWfOjof9n85suC0aXnGRA9dzAqLWuL6WcleQ7ScYn+Wyt9bZSyruHHv9kkm8leXWSnyV5IsnpvW4X6G7G1Mn5nV9/N+dO/KvMLA/lvuyTf/z1ezJjyoubLm1E9O2bcx+Ex6Znjki75qzOhTynHdh5zfrkSLWN37Dy9NrON6z00771q359zbZ0wdykP/YPgK5G4giM1Fq/lU6T4tnLPvms+zXJmSOxLegLDzzQud1vv523jWWLcsqKj2TSuM6bvFl5KKes+Eiy7EBv8tqs3z8I97O5C/LgD/42Z913Qj5xxlv74pSmJM18YByNjOxn/fwhv9sFc8f6vsH2kJMMqJG4iCewvVav7kw70+KLM6k+92J0k+pa34owFsxdkMx6SfL830ne+2NvyseQfrxobiPfsDIaGflsp1/bmfpFP38rjgvmQsdo5yS0xIgcgQG0kDd5MGr6+aK5Gx5dkXGb+S6xDY+u8FeQturn/HedIICB5r0H9KstvZnzJo+mfe41namP3HDusXndvJmZXNYlSSZPHJdT5s3si6/2rVMP2K7ltEA/5/9xF3SuC/RsrhMEMDA0MKBfeZMHo2b4V/tu6J+v9j3hwjxVJj1n2VNlUsafcGFDFbFV/Zz/cxckr12YjB8ak9MO7Mw71Q5gIDiFZKT10zm07DxlM8djjzQXg4RR1bdf7Tt3QT5/4115/YOXZs/8Ko9O3C9XT39HztiZWTIaGdnP+j3/5y5Ibrmic9/7LgaVnGRAaWBAE17wgtHZjjd5MGr69qt9k5zxnvOSz92QJJl++rU5Y2dvcLQysp/Jf+hvcpIB5RQSAAAAoPUcgQFNuP/+zu3++zdbx1jnr4rQn2QkQHdykgGlgQHJ6H8Qfuyxzq1fOgDDyUiA7uQkA8opJAAAAEDrOQKDbedwfQAAABqigQEAI0WjFwBgp3EKCTRh/PjOBINm2aI8vWJJNvzixjz9F4clyxY1XRFtJCMBupOTDChHYEAT/t2/a7oCGH3LFiXfODvjn17TmX9sRfKNszv35y5ori7aR0YCdCcnGVCOwABoo9Ov7bvTEe69+r8l65587sJ1T3aWAwDAVjgCA5pw772d2wMO2Pnb6rMPwYxdM8uq7VrOABvNjAQYi+QkA0oDA5rw6183XQGMujJtVrL6ns0vh2eTkWyN5jyDTk4yoDQwABgdx12QtV89K5Pq2o2L1pZJmXTcBQ0WxTbzgXHs8ZoB0GdcAwOA0TF3QSb9h08k4yd15qcd2Jl3AU8AALaBIzAAGD1zFyS3XNG576/DAABsBw0MaMLEiU1XANBeMhKgOznJgNLAgCYcfHDTFQC0l4wE6E5OMqBcAwMAAABoPQ0MaMI993QmAIaTkQDdyUkGlFNIoAlPPtl0BQDtJSMBupOTDChHYAAAAACtp4EBAAAAtJ4GBgAAANB6roEBTZg0qekKANpLRgJ0JycZUBoY0ITnP7/pCqA5p1/bdAW0nYwE6E5OMqCcQgIAAAC0Xk8NjFLKXqWU60opdw7dTt/MOgeWUr5XSrm9lHJbKeW/9LJN6Au/+EVnAmA4GQnQnZxkQPV6BMZ5SRbXWl+QZPHQ/KbWJzmn1npIkpcmObOUcmiP24Wxbe3azgTAcDISoDs5yYDqtYFxSpIrhu5fkeTUTVeotd5fa/3R0P3Hktye5IAetwsAAAAMkF4bGPvVWu9POo2KJDO6rVxKmZ3kxUl+0GWdd5VSlpRSlqxcubLH8gAAAIB+sNVvISmlXJ/keZt56APbs6FSyh5Jrk7yf9Vaf7Wl9WqtlyW5LEnmz59ft2cbAAAAQH/aagOj1nr8lh4rpTxQStm/1np/KWX/JA9uYb2J6TQvvlhr/bsdrhb6xa67Nl0BQHvJSIDu5CQDaqsNjK34epK3Jfnw0O01m65QSilJPpPk9lrrx3rcHvSHAw9sugKA9pKRAN3JSQZUr9fA+HCSE0opdyY5YWg+pZSZpZRvDa3z20nekuSVpZSlQ9Ore9wuAAAAMEB6OgKj1roqyXGbWX5fklcP3b8xSellO9B3fv7zzu3BBzdbB0AbyUiA7uQkA6rXU0iAHbFuXdMVALSXjAToTk4yoHo9hQQAAABgp9PAAAAAAFpPAwMAAABoPdfAgCbsvnvTFQC0l4wE6E5OMqA0MKAJBxzQdAUA7SUjAbqTkwwop5AAAAA73YO/WpMFn7opDz62pulSRlw/7xu0iQYGNOFf/7UzATCcjIS+tHDxnfnh8oez8Po7my5lxI36vslJBpRTSKAJTz/ddAUA7SUjoa+86PxvZ+36DRvnv/CDu/OFH9ydSRPG5aeXnNxgZb1rbN/kJAPKERgAAMBOc8O5x+Z182Zm8sTOR4/JE8fllHkzc8P7j224st71875BG2lgAAAAO82MqZMzZdKErF2/IZMmjMva9RsyZdKEzJgyuenSetbP+wZt5BQSAABgp3ro8bV58289P6cdfVC+dPPdWdlHF7s8+P5r86NJl2bP/CqP7rpfrr7/HUmOaLos6EsaGNCEKVOargCgvWQk9J1PHXlXcs1ZydK1uWTagclxFySZ33RZvVu2KGc88j+SPJkkmb7ugc78sjnJ3AU7b7tykgGlgQFN2H//pisAaC8ZCf1l2aLkG2cnT6/tzK++pzOf7NwP+aNh8cXJuiefu2zdk53lO3Pf5CQDyjUwAACAnafbh/yxbvWK7VsO9EQDA5pw552dCYDhZCT0l37+kD9t1vYtHylykgGlgQFNqLUzATCcjIT+0tSH/NFw3AXJxF2fu2zirkPX+NiJ5CQDSgMDAADYeZr6kD8a5i5IXrswmXZgktK5fe3CsX9tD2gpF/EEAAB2nmc+zC++uHPayLRZneZFv3zIn7ugf/YFWk4DAwAA2Ll8yAdGgAYGNGHatKYrAGgvGQnQnZxkQGlgQBP226/pCgDaS0YCdCcnGVAu4gkAAAC0ngYGNOGOOzoTAMPJSIDu5CQDSgMDAAAAaD0NDAAAAKD1NDAAAACA1tPAAAAAAFrP16hCE6ZPb7oCgPaSkQDdyUkGlAYGNGHffZuuAKC9ZCRAd3KSAdXTKSSllL1KKdeVUu4cut1iK7CUMr6U8s+llG/2sk3oCxs2dCYAhpORAN3JSQZUr9fAOC/J4lrrC5IsHprfkv+S5PYetwf94Wc/60wADCcjAbqTkwyoXhsYpyS5Yuj+FUlO3dxKpZRZSV6T5NM9bg8AAAAYQL02MPartd6fJEO3M7aw3v9Icm6SrR7nVEp5VyllSSllycqVK3ssDwAAAOgHW72IZynl+iTP28xDH9iWDZRSfi/Jg7XWW0opx2xt/VrrZUkuS5L58+fXbdkGAAAA0N+22sCotR6/pcdKKQ+UUvavtd5fStk/yYObWe23k7yulPLqJJOTTC2lfKHW+p93uGoAAABgoPR6CsnXk7xt6P7bklyz6Qq11v9Wa51Va52d5E1Jvqt5wcDbe+/OBMBwMhKgOznJgNrqERhb8eEki0opf5Dk7iRvSJJSyswkn661vrrH54f+5BcOwJbJSIDu5CQDqqcGRq11VZLjNrP8viTDmhe11u8n+X4v24S+sH5953ZCrz1EgD4kIwG6k5MMqF5PIQF2xF13dSYAhpORAN3JSQaUBgYAAADQehoYAAAAQOtpYAAAAACtp4EBAAAAtJ7L1kIT9t236QoA2ktGAnQnJxlQGhjQhOnTm64AoL1kJEB3cpIB5RQSaMJTT3UmAIaTkQDdyUkGlAYGNGH58s4EwHAyEqA7OcmA0sAAAAAAWk8DAwAAAGg9DQwAAACg9TQwAAAAgNbzNarQhP32a7oCgPaSkQDdyUkGlAYGNGHatKYrAGgvGQnQnZxkQDmFBJqwZk1nAmA4GQnQnZxkQGlgQBPuvrszATCcjAToTk4yoDQwAAAAgNbTwAAAAABaTwMDAAAAaD0NDAAAAKD1fI0qNOF5z2u6AoD2kpEA3clJBpQGBjRh6tSmKwBoLxkJ0J2cZEA5hQSa8MQTnQmA4WQkQHdykgGlgQFNWLGiMwEwnIwE6E5OMqA0MAAAAIDW08AAAAAAWk8DAwAAAGg9DQwAAACg9XyNKjRh5symKwBoLxkJ0J2cZEBpYEAT9tij6QoA2ktGAnQnJxlQPZ1CUkrZq5RyXSnlzqHb6VtYb89SylWllJ+UUm4vpbysl+3CmPf4450JgOFkJEB3cpIB1es1MM5LsrjW+oIki4fmN+cvk/zPWuv/keTIJLf3uF0Y2+67rzMBMJyMBOhOTjKgem1gnJLkiqH7VyQ5ddMVSilTk7w8yWeSpNb6VK310R63CwAAAAyQXhsY+9Va70+SodsZm1lnTpKVST5XSvnnUsqnSym7b+kJSynvKqUsKaUsWblyZY/lAQAAAP1gqw2MUsr1pZQfb2Y6ZRu3MSHJbyb561rri5P8Ols+1SS11stqrfNrrfP33XffbdwEAAAA0M+2+i0ktdbjt/RYKeWBUsr+tdb7Syn7J3lwM6utSLKi1vqDofmr0qWBAQAAALCpXr9G9etJ3pbkw0O312y6Qq31l6WUe0opL6q1/jTJcUn+d4/bhbFt1qymKwBoLxkJ0J2cZED12sD4cJJFpZQ/SHJ3kjckSSllZpJP11pfPbTeHyb5YilllyR3JTm9x+3C2Lbbbk1XANBeMhKgOznJgOqpgVFrXZXOERWbLr8vyaufNb80yfxetgV95Ve/6txOndpsHQBtJCMBupOTDKhej8AAdsQvf9m59UsHYDgZCdCdnGRA9fo1qgAAAAA7nQYGAAAA0HoaGAAAAEDraWAAAAAArecintCEgw5qugKA9pKRAN3JSQaUBgY0YfLkpisAaC8ZCdCdnGRAOYUEmrB6dWcCYDgZCdCdnGRAOQIDmvDAA53badOarQOgjWQkQHdykgHlCAwAAACg9TQwAAAAgNbTwAAAAABaTwMDAAAAaD0X8YQmzJ7ddAUA7SUjAbqTkwwoDQxowi67NF0BQHvJSIDu5CQDyikk0IRHHulMAAwnIwG6k5MMKEdgQBNWruzcTp/ebB0AbSQjAbqTkwwoR2AAAAAAraeBAQAAALSeBgYAALTBskXJxw9PLtqzc7tsUdMVAbSKa2AAAEDTli1KvnF2su7JzvzqezrzSTJ3QXN1AbSIBgY0Yc6cpisAaC8ZySBafPG/NS+ese7JznINDDYlJxlQGhjQhAn+6wFskYxkEK1esX3LGWxykgHlGhjQhFWrOhMAw8lIBtG0Wdu3nMEmJxlQGhjQBL90ALZMRjKAznn4lDxRd3nOsifqLjnn4VMaqohWk5MMKA0MAABo2Pvf98H87cz35d66TzbUknvrPrlq5vvy/nM/2HRpAK3h5CkAAGjYjKmTc8eMk3PR8sOyy/hxeerpDXnzjIPy1imTmy4NoDU0MAAAoAUeenxt3vxbz89pRx+UL918d1Y+tqbpkgBaRQMDAABa4FNvmb/x/iWnHt5gJQDtpIEBTfiN32i6AoD2kpEA3clJBlRPF/EspexVSrmulHLn0O30Laz33lLKbaWUH5dSvlxKcTIfg23cuM4EwHAyEqA7OcmA6nXUn5dkca31BUkWD80/RynlgCRnJ5lfaz08yfgkb+pxuzC2rVzZmQAYTkYCdCcnGVC9NjBOSXLF0P0rkpy6hfUmJNm1lDIhyW5J7utxuzC2PfJIZwJgOBkJ0J2cZED12sDYr9Z6f5IM3c7YdIVa671JPprk7iT3J1lda/1fW3rCUsq7SilLSilLVuoqAgAAANmGBkYp5fqha1dsOp2yLRsYui7GKUkOTjIzye6llP+8pfVrrZfVWufXWufvu+++27ofAAAAQB/b6reQ1FqP39JjpZQHSin711rvL6Xsn+TBzax2fJKf11pXDv3M3yX590m+sIM1AwAAAAOm11NIvp7kbUP335bkms2sc3eSl5ZSdiullCTHJbm9x+0CAAAAA6TUWnf8h0vZO8miJAel06h4Q6314VLKzCSfrrW+emi9DyV5Y5L1Sf45yRm11rXb8Pwrk/xihwtszj5JHmq6iJ3k4CQ/b7qInWS0X7fR3J4xOTb18+vWz/tmTI69bTWxvdFiPI5N/bxvxuTY3J4xOTaN5dft+bXWYdeU6KmBweaVUpbUWuc3XcfOUEr5da1196br2BlG+3Ubze0Zk2NTn79u/bxvxuQY21YT2xstxuPY1Of7ZkyOwe0Zk2NTP75uvZ5CAgAAALDTaWAAAAAAraeBsXNc1nQBO9HfNV3ATjTar9tobs+YHJv6+XXr530zJsfetprY3mgxHsemft43Y3Jsbs+YHJv67nVzDQwAAACg9RyBAQAAALSeBgYAAADQehoYAAAAQOtpYAAAAACtp4Exwkopny2lPFhK+XHTtfSilPK9UsoJQ/cvKaUsbLqmkdLP+7apfhmPSX+/bv28b5vqlzHZz69ZP+/b5hiT7dfP+7Y5xmT79fO+bapfxmPS369bP+/bpiY0XUAfujzJJ5J8vuE6enVhkotLKTOSvDjJ6xquZyT1875t6vL0x3hM+vt16+d929Tl6Y8x2c+vWT/v2+ZcHmOy7fp53zbn8hiTbdfP+7apy9Mf4zHp79etn/ftOXyN6k5QSpmd5Ju11sObrqUXpZS/T7JHkmNqrY+VUg5NclGSVUkWJ7k5nUB7KMkdtdYPN1Xr9tqGffvfz56vtV7VUKk965fxmBiTMSZbxXjsj/GYGJNjgTE5NhmT/TEm+2U8JsZk+mBMOgKDzSqlHJFk/yQP1VofG1p8cpL/p9Z6Qynl60keTXJtrfVTpZQx05Xdxn37+03mx+R/8H5iTBqTbWI8Go9tY0wak21jTBqTbWNM9seYdA0Mhiml7J/ki0lOSfLrUsqrhh76f5O8qZTy50n2TvLPQ/PfTfK9RordTtuxb5vO0yBj0phsE+PReGwbY9KYbBtj0phsG2Oyj8ZkrdU0wlOS2Ul+3HQdO1j7bkluSnLC0PzLk9y0yTrjk1yT5L8mefnQsquarn0k921L82NxGsvjcXtfN2NybExjeUwaj/03Hof2w5hs4WRMGpNtmwZ1TI7l8bi9r5sx2f7JNTB2gn46T+zZhvbrj5LsnuSv0znE6qJ0zhF7vNb6X5uqrVeb2bcVz56vtd7YXHW96dfxmBiTzVXXm34dk8bj2GVMjj3G5NhkTI5N/ToeE2Oyuep2nAbGCCulfDnJMUn2SfJAkgtrrZ9ptCgGlvFI2xiTtI0xSdsYk7SJ8UjbaGAAAAAArecingAAAEDraWAAAAAAraeBAQAAALSeBgYAAADQehoYAAAAQOtpYAAAAACtp4EBAAAAtJ4GBgAAANB6GhgAAABA62lgAAAAAK2ngQEAAAC0ngYGAAAA0HoaGAAAAEDraWAAAAAArTeh6QJor1LK1CR/n2SXJAcnuSPJmiT/vta6ocnaGEzGJG1jTNI2xiRtYjzSNsbk2FdqrU3XQMuVUo5O8oFa6ylN1wKJMUn7GJO0jTFJmxiPtI0xOXY5hYRtcXiS256ZKaXsXkq5opTyN6WUNzdYF4Nr0zE5p5TymVLKVQ3WxGCTk7SNnKRNZCRtIyPHKA0MtsWhSX78rPn/mOSqWus7k7yumZIYcM8Zk7XWu2qtf9BgPSAnaRs5SZvISNpGRo5RGhhsi5lJfvms+VlJ7hm6//TolwPDxiQ0TU7SNnKSNpGRtI2MHKM0MNgW30nymVLKK4bmV6TziycxhmjGpmMSmiYnaRs5SZvISNpGRo5RLuLJdiul7J7kE+lcsffGWusXGy6JAVdK2TvJf09yQpJP11r/tOGSGHBykraRk7SJjKRtZOTYoYEBAAAAtJ5DtgAAAIDW08AAAAAAWk8DAwAAAGg9DQwAAACg9TQwAAAAgNbTwAAAAABaTwMDAAAAaD0NDAAAAKD1NDAAAACA1tPAAAAAAFpPAwMAAABoPQ0MAAAAoPU0MAAAAIDW08AAAAAAWk8DAwAAAGg9DQwAAACg9TQwYBSUUsaXUv6ylHJbKeXWUsqcpmsCaBM5CdCdnAQNDBgt/y3JXbXWw5IsTPKehusBaBs5CdCdnGTgTWi6AOh3pZTdk/yHWutRQ4t+nuQ1DZYE0CpyEqA7OQkdGhiw8x2f5MBSytKh+b2SXN9cOQCtIycBupOTEKeQwGiYl+SCWuu8Wuu8JP8rydJSyu6llCtKKX9TSnlzoxUCNGteNp+Tc0opnymlXNVodQDNm5fN5+SpQ+8lrymlnNhohTAKNDBg55ue5IkkKaVMSHJikm8k+Y9Jrqq1vjPJ65orD6Bxm83JWutdtdY/aLQygHbYUk5+bei95NuTvLG58mB0aGDAzndHkpcO3X9vkmtrrT9PMivJPUPLn26iMICW2FJOAtCxtZw8P8mlo14VjDINDNj5vpzkN0spP0syN8n/PbR8RTpNjMT/RWCwbSknAejYbE6Wjo8k+Xat9UdNFgijodRam64BBtLQ1aQ/kWRNkhtrrV9suCSAViml7J3kvyc5Icmna61/2nBJAK1SSjk7yduS/DDJ0lrrJxsuCXYqDQwAAACg9Ry2DgAAALSeBgYAAADQehoYAAAAQOtpYAAAAACtp4EBAAAAtJ4GBgAAANB6GhgAAABA62lgAAAAAK33/wNlmeyogN4zIgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.errorbar(np.arange(n_periods * (len(het_inds) + 1) * n_treatments),\n",
    "             true_effect_params.flatten(), fmt='*', label='Ground Truth')\n",
    "plt.errorbar(np.arange(n_periods * (len(het_inds) + 1) * n_treatments),\n",
    "             param_hat.flatten(), yerr=((upper - param_hat).flatten(),\n",
    "                                        (param_hat - lower).flatten()), fmt='o', label='DynamicDML')\n",
    "add_vlines(n_periods, n_treatments, het_inds)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
